clc,clear
close all
%% STEP1:Set Parameters
driver.radius = 15e-3;
driver.diameter = 2*driver.radius;
driver.height = 30e-3;
driver.volume = pi*(driver.radius^2)*driver.height;
driver.magDirection = [0;0;1];
driver.position = [0;100e-3;0];
driver.attitude0 = [pi/2;0;0];

%% STEP2:Calulate B
rotB = [];

for theta = linspace(0,2*pi,360)
    attitude = [theta,0,0];
    rotB = [rotB [PMDip(Rot3(attitude),[0;-100e-3;0],driver);theta]];
end
theta = linspace(0,2*pi,360)
%% STEP3:Plot
fig = figure(1);
nImages = 360;
stepImages  = 10;
%设置曲线参数
h1 = animatedline;
h1.Color = 'r';
h1.LineWidth = 1.6;
h1.LineStyle = '-.';
h1.Marker = 'o';
h2 = animatedline;
h2.Color = 'b';
h2.LineWidth = 1.6;
h2.LineStyle = '-';
h2.Marker = '*'
%设置图窗参数
hold off%覆盖开关
% xlabel('Distance(mm)');
% ylabel('F(N)');
% xlabel('\theta(deg)')
% ylabel('\tau (N\cdot m)')
set(gca,'FontSize',20,'Fontname', 'Times New Roman');
axis equal
axis([min(rotB(2,:)),max(rotB(2,:)),min(rotB(3,:)),max(rotB(3,:))])

for idx = 1:stepImages:nImages
    plot([0,rotB(2,idx)],[0,-rotB(3,idx)],'r',[0,1e-3*cos(rotB(4,idx)+pi/2)],[0,1e-3*sin(rotB(4,idx)+pi/2)],'b')
    hold on
    
    axis([min(rotB(2,:)),max(rotB(2,:)),min(rotB(3,:)),max(rotB(3,:))])
    frame = getframe(fig);
    im{idx} = frame2im(frame);
    drawnow;
end
close;
%% 进行一个图片的保存
filename = ['我不想写程序了' '.gif'];
for idx = 1:stepImages:nImages
    [A,map] = rgb2ind(im{idx},256);
    if idx == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);
    end
end
function B = PMDip(rot3,pos,PM)
    xd = pos(1);
    yd = pos(2);
    zd = pos(3);
    Br = 1.465;
    V = PM.volume;
    nu0 = 4*pi*10^(-7);
    M = rot3*[0;0;1]*(Br*V)/nu0;
    r = [xd;yd;zd];
    D = norm(r);
    B = (nu0/(4*pi*D^3))*((3*r*r'/(D^2))-eye(3))*M;
end

function R = Rot3(attitude)
ax = attitude(1);
by = attitude(2);
cz = attitude(3);
Rz = [cos(cz) -sin(cz) 0;sin(cz) cos(cz) 0;0 0 1];
Ry = [cos(by) 0 sin(by);0 1 0;-sin(by) 0 cos(by)];
Rx = [1 0 0;0 cos(ax) -sin(ax);0 sin(ax) cos(ax)];
R = Rz*Ry*Rx;
end
